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Abstract: We study the QCD phase structure at high temperature and density adopting a histogram method. Because 
the quark determinant is complex at finite density, the Monte-Carlo method cannot be applied directly. We 
use a reweighting method and try to solve the problems which arise in the reweighting method, i.e. the sign 
problem and the overlap problem. We discuss the chemical potential dependence of the probability distribution 
function in the heavy quark mass region and examine the applicability of the approach in the light quark region. 
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1. Introduction 

It is important to study the QCD phase diagram including its quark mass dependence, as well as the temperature 
(T) and chemical potential (/i) dependences, for understanding the nature of the QCD phase transition at high 
temperature and density. The phase transition is expected to be first order when the up, down and strange quark 
masses are sufficiently large or small, and becomes crossover in the intermediate region between them. We expect 
a second order transition in the chiral limit of 2-flavor QCD and the nature of the chiral transition changes as the 
strange quark mass decreases. Moreover, the phase structure at nonzero [i is also expected to have a rich phase 
structure. To study the phase diagram, the probability distribution function (histogram) of a physical quantity 
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provides us with a useful method to determine the order of phase transitions. We expect that the distribution has 
two peaks at a first order transition when the volume is finite. Hence, we can identify the nature of the transition 
through the shape of the distribution function. In Sec. 2, we study the phase structure in the T, fi and quark 
mass parameter space measuring the probability distribution function when the quark is heavy, as the first step. 
We then extend this analysis to the light quark region and discuss the applicability of our approach in Sec. 3. 

2. Plaquette effective potential in the heavy quark region 

The probability distribution function for the plaquette value is defined by 

w(P;P,k,h)= fvU8(P{U]-P^ [detM(K,n)} Nt e m ' tt ^ P . (1) 

Here, we use the standard plaquette gauge action and N{ flavor Wilson quark action. P is the plaquette operator 
in the action and ft = 6/g 2 is the gauge coupling. M is the quark matrix and k is the hopping parameter 
corresponding to the inverse of the quark mass for a heavy quark. Applying the reweighting technique, we discuss 
the (ft, k, /^-dependence of the distribution function, 



w(P; /?,*,/*) = fVU6{P-P)[detM(K, li )] N <e 6pN «*° p = e(fi-Po) NMm P 
w(P; ft, ko,0) fT>U6(P-P)[detM(Ko,0)] N te 6 P° N »x° p 



det M(k,(i) 



-|JV f 



detM(K O ,0) 



(2) 



where {■ ■ -)p is the expectation value with fixed the plaquette P measured in a simulation at (kq,jj, = 0). We 
now define the effective potential as V e g(P ; ft, K, (a) = — \nw(P; ft, k, (/,). The (ft, K, /^-dependence is given by 



Veff (P; ft, k, ft) = V«ff (P; Po, /to, 0) - 6N sitc (p - ft ) P - In 



det M(k,/j.) 



-i N { 



detM(K ,0) 



(3) 



We study the heavy quark region where the hopping parameter n is sufficiently small, performing simulations at 
«o = 0. At the lowest order of the hopping parameter expansion, the quark determinants is evaluated as 



det M(k,,/i) 
det Af (0,0) 



= exp [288A" s it c K 4 P + 3 x 2 iYt+2 iV s 3 K JVt cosh (£\ |fi R + itanh (£\ flijl , (4) 



where Or, and Qi are the real and imaginary parts of the Polyakov loop. The first term proportional to P can be 
absorbed into the gauge action by a shift ft — > P + 48A r fK 4 . Because the term from the gauge action in Eq. (3) is 
linear, this term does not contribute to the appearance of the critical point at which the second derivative of V e s 
with respect to P vanishes. The complex phase of detM is induced by f2i, which causes the sign problem. 
We perform simulations of SU(3) pure gauge theory at five ft values in the range 5.68 - 5.70 on a 24 x 4 lattice. 
The details of the simulations are given in Ref. [1]. We combine the data obtained at five ft to enlarge the 
range of P. Measuring the histogram of P, i.e. the plaquette distribution function Eq. (1), we calculate the 
derivative of the effective potential. The /3-dependence of dV c g/dP at fixed k is controlled by dV c s/dP(ft) = 
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Figure 1. Derivative of V of f at nonzero K for iVj = 2. 
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Figure 2. Effevtive potential at f) c for each n. 



dVes/dP(0o) — 6A r s it c (/9 — j3o)- Since the additional reweighting term in this equation is constant, the shape of 
dV c ff/dP as a function of P does not change under a change of /3. Therefore, dV c g/dP is a good quantity to 
identify the order of phase transitions. If the transition is of first order, the effective potential has two minima. 
Then, dV c s/dP must be an S-shaped function and crosses the horizontal axis at three values of P. 
The re-dependence of the effective potential is investigated at fi = by the reweighting method up to the order re 4 . 
Using the data of w(P;fi,K = 0,/i = 0) and reweighting it, we evaluate dV c g/dP at nonzero re. Results near the 
critical re for JVf = 2 are plotted in Fig. 1. The S-shape structure becomes weaker as re increases and turns into a 
monotonically increasing function around re = 0.066. This behavior suggests that the first order phase transition 
at re = becomes a crossover at re « 0.066. We also show, in Fig. 2, the re-dependence of V e g at the transition 
point obtained by a numerical integration of dV c ff/dP. The double-well becomes shallower as re increases. 
We next consider the case of phase quenched finite density QCD, in which the complex phase term is neglacted 
in Eq. (4). In 2-flavor QCD, this corresponds to the case of the isospin chemical potential, /i u = —/id = \i- We 
find that, to the lowest order of the hopping parameter expansion, the effect of p is just to modify the hopping 
parameter as re — > re cosh 1 ' ip/T) in the theory at /i = 0. Therefore, the critical point at nonzero /i is given by 
k c P (m) — re cp (0)/cosh 1 ' iv *(^/r), where re cp (0) is the critical point at /i = 0. We plot the critical line in Fig. 3. 
The critical value decreases exponentially as fi/T increases. 

The complex phase of detM on the critical line is given by 8 = 3 x 2 Nt+2 N^ re^*(0) tanh(/x/T) fli. Because 
tanh(/i/T) < 1, the upper limit of the complex phase fluctuation can be estimated from the data of fli in the 
whole range of fi/T. It is found in Ref. [3] that the contribution from the complex phase to the location of the 
critical point is quite small in the analysis by quenched simulations on a 24 x 4 lattice. 

Adopting the phase quenched approximation, the application to JVf = 2 + 1 is straightforward because the 
difference from N{ = 2 is just the replacement of 2re JVt by 2re uc [ + ref*. We then find the critical (re u( j, re s ) satisfy 



2re^(M)cosh( Atud /T) + k?* (ji) cosh(/i s /T) = 2[k%= 2 (0)] Nt . 



(5) 



The critical lines in the re plane for up, down and strange are drown in Fig. 4 for p u d/T = 0-10 and fj, s /T = 0. 
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Figure 3. Critical point in the phase-quenched ap- 
proximation for iVf — 2. 



Figure 4. Critical line for each fi ud /T in N{ = 2 + 1 
QCD with fj, s /T = 0. 



3. Effective potential at fi ^ by phase quenched simulations 



Toward a study of the critical line in the light quark region, we propose a method based on the histogram 
method combined with the cumulant expansion of the complex phase of det M [2] . To handle the k and /j, 
dependence in a wide parameter range, we study the probability distribution function for the plaquette P and 
the quark determinant F(fj,) = N{ In | det M(^i)/det M(0)|. The combination of phase quenched simulations and 
the reweighting method is applied to evaluate the probability distribution function. Decomposing the quark 
determinant as (det M {fi)) Nf — e l6>(M) | detM(^)| JVf , the distribution function w(P,F; /3, fi) is calculated by 



/ 



w(P,F;P,v)= I VUS(P - P)S(F - F)e i9{f ' ) \det M(/x)| 



\N f 60N eite P 



«(c) 



(P,F) 



w (P,F,I3,h), (6) 



where wo(P, F; /3, fi) is the distribution function in the phase quenched simulation and {e' }(p,f) is the average of 
the phase factor with fixed P and F in the phase quenched simulation. Introducing the phase quenched effective 
potential Vb = — lnu>o, the effective potential V B s can be written as V a s — Vq — ln(e* }. 

In the study, there are two important problems. The first problem is that the effective potential in a wide range 
of (P, F) parameter plane is needed to investigate the shape of the effective potential. As is done in the previous 
section, we can enlarge the parameter range combining data obtained by different simulations. We perform 
simulations at several [i and combine these data using the reweighing method. Vo at /i is evaluated from the 
data at no by Vb(/x) = Vb(/io) — In R(/j,, fio), where R — wo(n)/wo(fio) = (| det M(fj,)/ det M(/io)|)(p i F)- Using this 
equation, we plot Vb at \±jT = 2.4 in Fig. 5. The data are taken with the RG-improved Iwasaki gauge action and 
Nf — 2 clover improved Wilson quark action on an 8 3 x 4 lattice. The data obtained at different simulation point 
are consistent. Combining these data, we evaluate Vb in a wide range of F. 

The next problem is the sign problem. If the sign of e lS changes frequently, the error of (e 1 }(p,f) becomes larger 
than the mean value. To aboid this problem, we calculate {e l0 )(p,F) by a cumulant expansion, (e ! )(p,f) — 

exp[i(6) c - {0 2 ) c /2 - i((9 3 ) c /3! + (6> 4 ) c /4! H ], where (8 n ) c is the n-th order cumulant of 9. In this equation, the 

odd order cumulants cause the sign problem, but these terms must be zero from the symmetry under /i to — /j,. 
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Figure 5. Effective potential V at (i/T = 2.0 eval- 
uated at 3 simulation points. 
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Figure 6. Distribution of the complex phase of the 
quark determinant. 



Once we remove the odd terms, the sign problem changes to the convergence problem of the cumulant expansion. 
When 9 is small, this expansion converges. However, if the distribution of 8 is almost uniform in — n < 6 < n, 
the expansion dose not converge. We adopt an alternative definition of the phase in the range — oo < 6 < oo. We 
measure the derivatives of lndet M with respect to fi, and define the phase by integrating the derivatives over /i, 



8(v) = TVflm [lndet MO)] = N f 



H/T 



Im 



9(lndetM(/x)) 



(7) 



We plot the distribution of 9 in Fig. 6. The distribution is well approximated by a Gaussian function. The 
dashed line is the fitted Gaussian function. For such a case, the cumulant expansion converges even when the 
fluctuation of 8 is large, since (e* } is given only by the second order cumulant if the distribution is Gaussian 
exactly. Therefore, the result of the 9 distribution suggests that the estimation of the complex phase factor by 
the cumulant expansion works well even in the high density region where the complex phase fluctuation is rather 
large. The details of the results shown in this section are given in Ref. [2]. 



4. Summary 



We studied the QCD phase diagram in the heavy quark region. The critical surface which separates the first 
order region and the crossover region in the (K u d, k s ,/i) space is calculated performing simulations in the heavy 
quark limit. We moreover discussed a new method to study the critical surface in the light quark region. 
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